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ABSTRACT 



Aims. The main aim of the paper was performing test of our (chemical and kinetic) codes, which will be used during self-consistent 
modelling of dynamics and chemistry in the winds from C-rich AGB stars. 

Methods. We use the thermodynamical equilibrium code to test the different databases of dissociation constants. We also calculate 
the equilibrium content of the gas using the kinetic code, which includes the chemical network of neutral-neutral reactions studied by 
Willacy & Cherchneff (1998). The influence of reaction rates updated using the UMIST database for Astrochemistry 2005 (UDFA05), 
was tested. 

Results. The local thermodynamical equilibrium calculations have shown that the NIST database reproduces fairly well equilibrium 
concentrations of Willacy & Cherchneff (1998), while agreement in case of Tsuji (1973) dissociation constants is much worse. 
The most important finding is that the steady state solution obtained with the kinetic code for reaction network of Willacy & 
Cherchneff (1998) is different from the thermodynamical equilibrium solution. In particular, CN and C2, which are important 
opacity sources are underabundant relative to thermodynamical equilibrium, while O-bearing molecules (like SiO, H2O, and OH) 
are overabundant. After updating the reaction rates by data from the UDFA05 database consistency in O-bearing species becomes 
much better, however the disagreement in C-bearing species is still present. 
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1. Introduction 

Asymptotic Giant Branch (AGB) stars develop a molecular 
dusty wind, driven by radiation pressure on dust grains 
condensating in the warm molecular layers. The new formed 
molecules and dust are deposited then into interstellar medium 
(ISM) and contribute to the global matter budget of galaxies. 
Therefore, the determination of molecular and dust contents in 
the AGB winds is important for understanding of their further 
re-processing in the ISM before, eventually, matter ends up in a 
new star forming regions. 

Chemical evolution of gas around AGB stars takes 
place at such different conditions like thermodynamical 
equilibrium in stellar photospheres, non-equilibrium (including 
possible shock-induced) chemistry, dust formation and photon 
dominated chemistry in circumstellar envelopes. 

Observations of circumstellar envelopes at infrared and 
millimetre range demonstrate that they are very rich in 
molecules. In the case of the well-studied C-rich AGB star 
IRC+10216 the number of identified molecules exceeded 60 
jOIofssoril I2OO6). The majority of the observed molecules 
are carbon-chain molecules such as the cyanopolyynes, the 
hydrocarbon radicals and carbenes, as well as organo-sulphur 
and organo-silicon molecules. Interferometric and single dish 
observations show that carbon-chain molecules are distributed 
in the hollow shells around central star Such distribution 
has been interpreted as arising from photochemistry of parent 
molecules (e.g. C2H2 and HCN) flowing away from the star (see 



e.g. reviews by iGIassgoIdl [19961 'Millar"2003'). The new space 
and ground-based projects (e.g. Stratospheric Observatory for 
Infrared Astronomy - SOFIA, The Atacama Large Milimeter 
Array - ALMA or Herschel Space Observatory) will allow us to 
investigate the circumstellar molecular inventory in much more 
details by using, not yet fully exploited, sub-millimetre range. 

With the aim to prepare tools for analysing the new sub- 
millimetre observations of AGB stars we are in a process 
of developing a self-consistent treatment of dynamics and 
chemistry in the circumstellar envelopes. The first effort to 
build the chemical network of reactions for conditions close to 
thermodynamical equilibrium (hi gh density and tempera t ure) in 
C-rich environment was done bv lWillacv & Cherchneff! (Il998h 
(hereafter WC98). Since the considered densities are high and 
effective temperature of the star is low, only neutral-neutral 
bimolecular and termolecular reactions were considered. WC98 
investigated chemistry of sulphur and silicon in the inner wind 
of IRC -1-1 02 16. Their non-equilibrium calculations started from 
chemical composition derived at the local thermodynamical 
equilibrium (LTE) conditions. It was found that the parent 
molecules C2H2 and CO are unaffected by the shocks, but 
the abundances of daughter molecules can be significantly 
altered by the shock-driven chemistry. Their model indicated 
that molecule production by the shock-driven chemistry close 
to the photosphere can be significantly more important than 
the LTE chemistry. Recently, Cherchneff (2006) investigated the 
chemistry and composition of the quasi-static molecular layers 



of AGB stars using the same semi-analytical model of shock 
dynamics as in WC98. 

In modelling of stellar atmospheres it is commonly assumed 
that conditions of local thermodynamical equilibrium prevails. 
Molecular content is then defined by dissociation constants. 
In the outflow, in the presence of shocks, kinetic equations 
must be solved for some chemical network. When local 
thermodynamical equilibrium is valid, at some position in 
the stellar atmosphere, both approaches should produce the 
same molecular concentrations. If not, then an inconsistency 
will appear, darkening the interpretation of obtained results. 
The aim of our paper is to clarify this point for future 
application of available chemical networks in the non- 
equilibrium computations of chemistry in the inner layers of 
carbon-rich outflows. 

In this paper we present our steady-state chemical code 
and test (some of) the available chemical databases using 
for comparison work of WC98. In the next paper (Schmidt 
et al. in preparation) we will discuss our approach to self- 
consistent treatment of dynamics and chemistry. The first 
prel iminary results of such computations were presented by 
iPule cka et al. (2005a) (using the TITAN hydrodynamical code) 
and IPulecka et al.i (i2005h) (using the FLASH hydrodynamical 
code). 



2. Test of chemical codes: The thermodynamical 
equilibrium 

In the stellar photospheres the dynamical timescales are much 
longe r than the timescales of molecular formation ( iGlassgoidI 
119961) . Therefore, to determine local molecular contents it 
is suflicient to apply the thermodynamical equilibrium (TE) 
conditions. To compute molecular composition in the local TE 
conditions one needs to know the temperature (T), the total gas 
number density («), or equivalently the total gas pressure (p), 
and the initial elemental abundances. We have constructed such 
TE code and present results of the performed tests. 



2. 1 . The thermodynamical equilibrium code 

The equilibriu m chemistry code is based on the Russell's 
jRusselll Il934l) approach. This method uses the fictitious 
pressure Pf (X), which is the pressure exerted by element X if all 
the gas were in the neutral form of the atomic species X. Hence 
the fictitious pressure of element X is given by: 



;:'f(X) = px -H ^ Q'x.yPy, 



(1) 



where: py - partial pressur43 of molecular species Y, which 
contains element X, q-x.y - the stoichiometric coefficient 
indicating the number of element X involved in molecule Y. For 
example, in the case of hydrogen we have: 



Pf (H) = ph+ 2/?h2 + PCH + 2pcH2 + 



(2) 



(plus all examined hydrogen-bearing molecules). 

Partial pressure of compound X is related to its concentration 
«x [cm^^] by the ideal gas law px - iixkT [dyn cm"^], where 
k - Boltzmann constant [erg/K], and T - temperature [K]. Using 



' The partial pressure of species Y is the pressure of the gas if no 
other species were presented in the medium. 



the partial pressures, the dissociation constant (Kp) of molecule 
AB (AB ^ A -I- B) can be expressed as: 



/^p(AB) 



PaPb nA kT tiB kT riA hb 



Pab 



«AB kT 



kT, 



(3) 



where: p^, ps, Pab - partial pressures of atoms A, B, and 
molecule AB, ha, n^, hab - their concentrations. Therefore by 
employing the Eq. [3] the fictitious pressure of hydrogen can be 
rewritten to: 



PfiH) ^ph+2 



Pi 



PhPc 



+ 2- 



pipe 



K^{W2) ^p(CH) K^{CW2) 



(4) 



Note, that dissociation constants can be determined from the 
differences of the corresponding Gibbs free energie^ Since, the 
fictitious pressure reflects the conservation of mass for element 
X, then /?f(X) for any other atomic species X can be expressed 
in terms of the fictitious pressure of hydrogen /?f (H): 



Pf(X) = A(X)pf(H), 



(6) 



where: A(X) = — is abundance of element X relative to 
hydrogen. 

To close this set of nonlinear equations for fictitious 
pressures we need only the equation for the total gas pressure, 
which is obtained by summing the partial pressures of all 
examined species. This final set of equations is solved using 
the Netwon-Raphson's method. The solution gives the partial 
pressure of atoms, which are used then to calculate molecular 
partial pressures (i.e. in fact their concentrations) by employing 
Eq.E] 



2.2. The equilibrium case of WC98 

In a first test of our thermodynamical equilibrium code we aimed 
to reproduce the thermodynamical equilibrium concentrations 
of species given in Table 3 of WC98. This table contains 
equilibrium abundances determined at distance of 1.2 R» and 
physical conditions specified in the entry of their Table 2 
(J = 2062 [K] and n = 3.68x10^4 [^m ^] what corresponds 
to /7= 1.05x10^ [dyn cm"^]). To reproduce the equilibrium 
molecular content at these conditions we need to know the initial 
elemental abundances and the dissociation constants. 

The initial elemental abundances were determined by 
summing abundances of all species (see Table 3 of WC98) 
containing a given element. The obtained values are given in 
Col. 2 of Table [T] where abundance of element X is given 
by e(X) - log'-j^ + 12, and nx is concentration of element 
X. For compa rison, the solar abundance s of H, He, C, O, N, 
Si, and S from'Grevesse & Sauvall (Il998h (hereafter GS98) and 
lAsplund et al. (2005) (hereafter A05) are shown in Cols. 3 
and 4, respectively. The initial element abundances assumed by 
WC98 are essentially solar ones except of C/O ratio, which was 
assumed to be about 1.5. 



- Here the sign ° marks conditions when temperature is 273.15 K 
and the absolute pressure is 1 atm (i.e. 1.01325 bar, 1.01325X 10* [dyn 
cm-^], 101.325 [kPa]). (AG°) of products (A, B) and reactant (AB) 
according to: 



Kp{AB) = C^" exp\- 



AG° 
RT 



(5) 



where: R = 8.32441 x 10' [erg K"' mor'] is the ideal gas constant, C = 
1.01325X 10* [dyn cm"^] and Ay is the stoichiometric factor of a given 
dissociation reaction. 
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Fig. 1. The ratio of NIST and WC98 equilibrium concentrations. The order of species (except of elements, which are placed at the 
right side of the figure) follows decreasing concentrations of WC98. 



Table 1. The initial abundances of elements: e(X) = log— + 12. 



element 


WC98 


GS98 


A05 


(1) 


(2) 


(3) 


(4) 


H 


12.00 


12.00 


12.00 


He 


10.99 


10.93 


10.93 


C 


8.98 


8.56 


8.39 


N 


8.06 


8.05 


7.78 


O 


8.79 


8.93 


8.66 


Si 


7.49 


7.55 


7.51 


S 


7.19 


7.33 


7.14 



The dissociation constants were obtained using species 
concentrations given in Table 3 of WC98 by employing Eq. |3] 
Note that the dissociation constants derived in such a way are 
valid only for the conditions (temperature) specified above. As 
expected, with such ^'s and initial elemental abundances given 
in column 1 of Table[l] the obtained equilibrium concentrations 
are in a very good agreement with values given in Table 3 of 
WC98. However, after discussion with K. Willacy, we realized 
that there are some misprints for the concentrations of C2 and 
SiH4 in Table 3 of WC98. Therefore, during further computation 
we have implemented the original AG°'s (kindly provided by K. 
Willacy) from work of WC98. These AG°'s are collected now in 
Table |2] (Cols. 3 and 8), and were used to determine equilibrium 



species concentrations which will be referred hereafter as WC98 
concentrations. The molecules listed in this table (with ordinal 
numbers in Cols. 1 and 6) are organised in order of decreasing 
WC98 concentrations. 

2.3. Comparison of TE concentrations based on WC98, 
NIST and Tsuji databases 

As was described in Sect. 12.21 determination of 
TE concentrations at specified conditions requires knowledge 
of KpS (or equivalently AG°'s). The most comprehensive 
database, which is publicly available and allows to compute 
dissociation constants is that collected by the National Institute 
of Standards and Technology (NIST). The NIST chemistry 
WebBool|3 provides thermoc hemic al data, like enthalpy of 
formation {H°) and molar entropy (S °) for different temperatures 
at standard conditions. These quantities allow determination of 
the Gibbs free energy (G°) for all examined species according 
to: 

G°^H°-TS°. (7) 

Thus, we can determine the difference of Gibbs free energies 
(AG°) for products and reactants, of the direct dissociation of 

^ http://webbook.nist.gov/chemistry/ 
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Fig. 2. The differences of AG° between database (DB: DB = NIST or Tsuji) and WC98. Results are marked by circles for the NIST 
database and by pluses for the Tsuji database. The order of species is same as in Fig.[T] 



molecule into single atoms, and use them to find the dissociation 
constants from Eq. |5] 

We computed concentrations of species according to the 
method described in Sect. |2.1| using the dissociation constants 
obtained from the NIST database for T = 2062 [K]. The 
NIST database does not contain thermochemical data for some 
molecules: C3H - #9 (number appearing after symbol # is the 
ordinal number of a given molecule used in Fig.lTJ, C4H2 - #16, 
HCS - #25, C4H - #26, HCSi - #27, HNSi - #28, SiHj - #29, 
SiCHa - #30, C3H2 - #31, SiHj - #37, so we have used AG°'s 
from WC98, instead. The comparison between NIST and WC98 
concentrations is shown in Fig.[T] The order of species (except 
of elements, which are placed at the right side of the figure) 
follows decreasing concentrations of WC98. The consistency is 
fairly good (almost all results agree within about 7 %, with the 
exception of SiC2 - #11 concentration, for which difference is 
about 22%. 

From chemical point of view a better illustration of 
discrepancies between different databases of thermochemical 
data is comparison between AG°'s. Therefore, in Cols. 4 and 
9 of Table|2] we provide the AG°'s derived from the NIST 
database. We left empty spaces for species with missing data 
in that database. In Fig.|2] we plotted the differences between 
AG°' from NIST and WC98 (circles) for each species. Again, 
the order of species (except of elements, which are placed at 



the right side of the figure) follows decreasing concentrations of 
WC98. For molecules without data in the NIST database there 
are no corresponding symbols in Fig.|2] As one can see energies 
agrees within about 2 [kJ mol"'] with exception of SiC2 for 
which discrepancy is about 5 [kJ mol"']. Another widely used 
database, especially in the modelling of stell ar atmospher es, is 
the set of dissociation constants compiled bv lTsuiil (Il973[) . The 
thermal dependence of /Tp's in this database is given by: 

logK^^cct + cie + ciO^ + c^e^ + CAei^, (8) 

here: Ci i=o.i 4 - are coefficients derived by Tsuji and - ^y^. 

This equation allow us to compute dissociation constants at the 
temperature of interest (T - 2062 K). The corresponding AG°'s, 
computed from Eq.|5] are collected in Cols. 5 and 10 of Table |2] 
and compared with WC98 and NIST values in Fig.|2] (pluses). 
The Tsuji's database does not contain seven molecules among 
46 considered: C4H2 - #16, HCS - # 25, C4H - #26, HCSi - 
#27, HNSi - #28, SiCH2 - #30, and , C3H2 - #31 (empty spaces 
in Table |2] and in Fig.|2]i. In general, agreement between AG°'s 
from Tsuji and WC98 is worse than that in case of NIST. The 
largest differences are seen for: CS - #7, CN - #10, SiC2 - #1 1, 
SiH - #12, C3 - #13, CH3 - #15, CH2 - #17, SiN - #23, NH2 - 
#35, HCO - #36, and NS - #39. AG°'s for other molecules differ 
less than about 4 [kJ/mol]. 



4 



Table 2. AG° in [kJ/mol] for T = 2062 K and = 1 atm in case of WC98, NIST, and lfsuiil ( [19731) . 



# 


species 


WC98 


NIST 


TSUJI 


# 


species 


WC98 


NIST 


TSUJI 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


1 


H2 


207.6 


206. 1 


205.6 


24 


SIC 


215.4 


214.5 


213.9 


2 


CO 


801.6 


801.6 


800.6 


25 


HCS 


463.0 






3 


C2H2 


878.4 


878.0 


874.4 


26 


C4H 


1279.8 






4 


N2 


684.3 


684.4 


685.1 


27 


HCSi 


386.0 






5 


C2H 


676.0 


675.8 


673.0 


28 


HNSi 


476.7 






6 


HCN 


770.0 


769.8 


773.5 


29 


SiH2 


150.4 




150.4 


7 


CS 


453.8 


453.2 


501.1 


30 


SiCH2 


555.8 






8 


SiS 


372.7 


371.4 


370.8 


31 


C3H2 


1038.2 






9 


C3H 


1049.7 




1050.0 


32 


NH 


109.9 


109.9 


109.3 


10 


CN 


507.0 


507.1 


519.1 


33 


H2O 


455.4 


455.4 


455.4 


11 


SiC2 


740.9 


736.2 


733.6 


34 


NH3 


464.7 


464.6 


464.7 




c;u 
3irl 


OA C 




103.0 




NH2 


277. 1 


277.2 


inn fi. 


13 


C3 


794.3 


793.5 


851.1 


36 


HCO 


662.7 


662.5 


716.0 


14 


HS 


152.5 


152.0 


151.1 


37 


SiH3 


280.9 




281.0 


15 


CH3 


516.5 


516.2 


524.9 


38 


S2 


183.7 


183.4 


186.9 


16 


C4H2 


1550.3 






39 


NS 


262.4 


261.7 


282.0 


17 


CH2 


318.8 


318.6 


282.5 


40 


OH 


212.1 


212.1 


212.6 


18 


CH 


131.7 


131.6 


129.1 


41 


H2CO 


763.3 


763.4 


762.3 


19 


SiO 


537.6 


536.9 


538.0 


42 


CO2 


1024.3 


1024.1 


1023.1 


20 


C2 


346.4 


346.2 


343.9 


43 


SiH4 


335.1 


333.9 


336.9 


21 


H2S 


286.8 


286.2 


288.4 


44 


NO 


395.5 


395.3 


395.7 


22 


CH4 


(,61.1 


667.4 


666.6 


45 


SO 


279.3 


278.7 


279.2 


23 


SiN 


317.3 


316.6 


261.4 


46 


0. 


234.8 


234.8 


235.0 
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Fig. 3. The ratio between Tsuji and NIST equilibrium 
concentrations at conditions specified in Sect. 12. 21 - black circles. 
Red circles mark the ratio, when AG° for CS in the Tsuji 
database was replaced by the value from WC98. The order of 
species is same as in Fig.[T] 

The influence of differences between the Tsuji and NIST 
AG°'s on species concentration can be seen in Fig.[^ - black 
circles. The order of species is same as in Fig.[T] We plotted 
there the logarithm of Tsuji concentrations relative to those of 
NIST (using the WC98 values of AG° in case of molecules with 
missing data). The consistency is rather poor. Concentrations 
for all molecules with considerable discrepancies in AG°'s 
are significantly different, but additionally some other species 
reached well different equilibrium state. For example, for all 
sulfur-bearing species, except of CS - #7 (i.e. SiS - #8, HS 
- #14, H2S - #21, HCS - #25, S2 - #38, NS - #39, SO - 
#45, S - #53), we see very large discrepancies (even by factor 
of 84! for S2), in spite of the fact that differences between 
AG°'s are significant only in case of CS and NS. Among other 
molecules with significant discrepancy in concentrations are 

Figs. 3-5 and 7 are available only in the online version. 



Fig. 4. The ratio between Tsuji (with AG° for CS replaced 
by the value from NIST) and NIST equilibrium concentrations 
at conditions specified in Sect. 12. 21 - black circles. Red circles 
mark the ratio, with additional replacement of AG° for CS, CN, 
SiC2, SiH, C3, CH3, CH2, SiN, NH2, HCO, and NS in the Tsuji 
database by the values from WC98. The order of species is same 
as in Fig.[T] 



those containing C: C3 - #13, CH2 - #17 and HCO - #36; two 
silicon-bearing molecules: SiN - #23 and SiH - #12; and radical 
NH2-#35. 

2.4. The influence of changes in the AG° 's on the TE 
concentrations 

Here we test the influence of changes in AG° on equilibrium 
molecular content of examined gas. 

In the first test we only replaced the Tsuji AG° for 
CS (this molecule has the largest equilibrium concentration 
among S-bearing molecules) by the smaller value from the 
NIST database. The correlation between the Tsuji and NIST 

concentrations becomes much better (see red circles in Fig.[J|), 
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Fig. 6. Comparison between species concentrations in [cm ] derived by our kinetic code for case 1 - crosses, and for case 2 
- circles, at conditions specified in the 1" entry of Table 2 in WC98. Squares mark concentration for case 2, with rate constants 
updated using the UDFA05 database. The order of species is same as in Fig.[T] 
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Fig. 5. The ratio between Tsuji and NIST equilibrium 
concentrations at conditions specified in Sect. 12. 21 - black circles. 
Red circles mark the ratio, when AG° for C2H2 in the Tsuji 
database was replaced by the value from NIST. The order of 
species is same as in Fig.[T] 



especially for family of all S-bearing molecules. Decrease in 
concentration of CS - #7 releases atomic carbon and sulfur, 
and the excess of sulfur increases concentrations of all S- 
bearing species. However at the same time, the increase in 
SiS - #8 concentration cause the decrease of concentration 
for all other silicon-bearing species, even for carbon-bearing 
ones: SiCi - #1 1 and SiC - #24. Therefore, it seems that SiS 
molecule controls the silicon chemistry in the examined case. 
The considerable increase in concentration of sulfur monoxide 
(SO - #45) is caused by both: an access to sulfur due to decrease 
of ties and an access to oxygen due to decrease of nsio (SiO - 
#19). Furthermore, enrichment of the gas in carbon results in 
increasing concentrations of all molecules composed of only C 
or carbon and hydrogen. Note, that decreasing of AG° for NS - 
#39 (molecule with the second largest difference between AG° 
in the NIST and Tsuji databases) would result in much better 
agreement between corresponding concentrations only for this 
molecule). 

As a starting point for the second test we chosen 

concentrations marked by the red circles in Fig.[^. They are 

repeated now by black circles in Fig.|^. The order of species 

is same as in Fig.[T] Red circles in Fig.|^ show now the 
Tsuji concentrations relative to the NIST ones obtained after 
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additional replacement of AG°'s by the NIST values for all 
molecules, which show significant discrepancy in the Gibb's 
energies (i.e. for CN - #10, SiCa - #1 1, SiH - #12, Q - #13, CH^ 
- #15, CHz - #17, SiN - #23, NHj - #35, HCO - #36, and NS - 
#39 - see Fig. |2] and Table|2]). As we can see, changes in AG°'s 
for these molecules alter almost only concentrations for these 
species. The exception is C3, where atomic carbon released by 



its decreasing abundance (due to the smaller AG° in the NIST 
database) is bound also into other carbon-bearing molecules: 
C3H - #9, SiC2 -#11, C4H2 - #16, C4H - #26, and C3H2 - #3 1 . 

Finally, we present results, which show how relatively small 
changes of AG° for abundant molecule like acetylene (C2H2 - 
#3) alter the equilibrium abundances of other molecules. AG° 
for acetylene from the Tsuji database is only about 3.6 [kJ/mol] 
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smaller than that from NIST (see Table|2]i. In Fig. |^ we present 
the ratio between the Tsuji and NIST equilibrium concentrations 
by black circles, while red circles show these ratios computed 
with AG° replaced only for acetylene by the larger value from 
the NIST. As we can see, even so small change in AG° for this 
parent and abundant molecule alter significantly the obtained 
concentrations for almost all other species. With larger AG° for 
acetylene more carbon is used for formation of this molecule 
and therefore concentration of other C-bearing molecules are 
significantly decreased. This affects also concentrations of other 
species. In case of silicon-bearing molecules all silicon released 
from SiC2 is bounded into SiS and SiO. The less sensitive to the 
change of AG° for acetylene turned out to be such molecules as: 
SiH - #12, HCS - #25, SiHz - #29, HCO - #36, SiH^ - #37, 
H2CO - #41, and SiH4 - #43. 

3. Test of chemical codes: The chemical kinetics 

In the outer part of the stellar atmospheres one can expect that 
non-TE effects are becoming more and more important and non- 
TE methods to determine properly concentrations of chemical 
species are necessary. We have developed a kinetic code, 
which can treat time-dependent chemistry and describe here the 
performed tests. Now, however, instead of using dissociation 
constants (or equivalently AG°'s) we need to consider an 
appropriate reaction network. In this work we consider reaction 
network of WC98. 

3.1. The kinetic code and reaction network 

To find molecular concentrations we have to consider all 
chemical reactions from the given network, which lead to 
formation and destruction of all molecules under study. The 
evolution of concentration with time for species X can be written 
as: 



dt 



(9) 



where: the Fx stands for reactions responsible for production 
of species X, and Dx for its destruction. In the developed 
code the rate equations, defined for each considered reaction, 
are accumulated. The rate equation describes the rate of 
disappearance of each of the reactants and on the other hand 
the appearance rate of each of the products. For example, for 
reaction A + B C + D the rate equation is defined as: 



driA 
dt 



driB 
dt 



dnc 
IT 



driu 

~dr 



— k iPj^ /ig, 



(10) 



where: «a, hb - are the concentrations of reactants, rtc, «d - 
are the concentrations of products, and the parameters y and 
6 characterise order of reaction with respect to each reactant. 
These parameters are derived experimentally and their sum 
determines the reaction order. The coefficient k in the above 
equation is known as a reaction rate constant and determines 
reaction dependence on the temperature. In our model the order 
of considered reactions was always the same as the reactants 
number, i.e. 7 and 6 were equal 1. 

We search for the steady state solution of the above set of 
rate equations using DVODE solver for systems of ordinary 
differential equation written at Lawrence Livermore National 
LaboratorjQ DVODE solves the stiff system of differential 



equations of first order, just like systems which are often met 
in the chemical kinetics. The solver must be updated with 
the subroutine providing the net reaction rates and eventually 
subroutine providing Jacobian of the system. These two routines 
are constructed automatically by searching the network for given 
molecules and preparing set of corresponding rate equations for 
given range of temperature. The WC98 reaction network used 
in our study includes 53 species (7 elements and 46 molecules 
- given in Table[T]and Table|2] respectively). The molecules are 
composed of six elements: H, C, N, O, Si, and S. In Table 5 of 
WC98 there are data for 352 reactions. Among them there are 
139 reactions in both directions (hereafter "two-way") and 213 
reactions which do not have the backward reaction documented 
(hereafter "one-way", marked by " <->" in their table). However 
there are some miss printings in this table: 

1. The "two-way" reaction number 53 (H2 -H N2 ^ NH + NH) 
has rate constant documented as being zero. We computed 
this missing rate constant (by method described below - 
see Eqs.[T2] and \T3[ using equilibrium constant and the 
rate constant for its backward reaction (i.e. for reaction 
number 269: NH + NH -> H2 + N2). Note, however, that 
such approach causes that reaction 269 must be treated as 
"one-way" process. After this change number of "one-way" 
reactions in network increased by one to 214, while number 
of "two-way" reactions decreased to 137; 

2. The reaction number 78 (C + NS -> CS -H N) is indicated 
as a "two-way" reaction (there is no " <->" in the table), but 
this reaction does not have backward reaction documented 
and, in fact, should be traeated as "one-way" process. This 
change caused that number of "one-way" reactions increased 
to 215, while the number of "two-way" reactions decreased 
to 136; 

3. The reaction number 217 (S-t-N2 ^ S + NS) is incorrect. 
It seems that the correct version of this reaction is: Sh-N2 
— > N + NS. However, such reaction is already present in 
the network, therefore we simply deleted reaction number 
217 from the list. This way, number of "one-way" reactions 
decreased by one to 214; 

4. The coorect version of reaction 279 (SiH + NO SiO + H) 
should be SiH -h NO ^ SiO + NH. This does not change the 
statistics of rections. 

In summary, we have 564 reactions in the network according 
to Table 5 of WC98: 214 in "one-way" + 214 backward - 
determined as described below +136 "two-way" reactions. 

The rate constants for forward and backward reactions 
documented in Table 5 of WC98 were computed via equation 
in Arrhenius form: 



^(300) 



(-1)' 



(11) 



where: kf and kh are the rate constants of forward and backward 
reactions, respectively. A, /3, and are the Arrhenius parameters 
given in that table. In case of missing backward reactions, 
WC98 calculated the backward rate constant ki, from the 
thermodynamics of the reaction (see Eq. 4 in WC98). In our 
computations the backward rate constant k^ for specific reaction 
and given temperature was derived from the ratio of rate constant 
for forward reaction, and the reaction equilibrium constant Ki, 
i.e.: 



kf/K,. 



(12) 



http : ://w w w.Unl . gov/CAS C/odepack/ 



On the other hand, the equilibrium constant for given 
reaction can be obtained from the ratio of equilibrium 
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concentrations for products and reactants. For example, for 
reaction AB + C — > AC + B the corresponding equilibrium 
constant is given by: 



»AC«B 

"Afinc ' 



(13) 



where: «ab, "c, "ac, "b - equilibrium concentrations of 
reactants and products, respectively. Therefore after inserting Kr 
given by Eq. ( fT3] l into Eq. (fTZt we have: 



"ABnc 

f ■ 

"acWb 



(14) 



Thus, by computing missing backward rate constants, we have 
modified the original network of WC98 and use it to test our 
kinetic code. 



3.2. Comparison of TE concentrations derived witli WC98 
and UDFA05 reaction networks 

The chemical timescales in the dense and warm part of the 
envelope are very short, so the concentrations of species should 
reach the steady state values very quickly. Therefore, to test 
our kinetic code we decided to reproduce the equilibrium 
concentrations of the 53 species discussed above. In such case, 
as an input to the code we need the initial concentrations of 
elements only (values are taken from Table[T]), the total gas 
density and its temperature (again we have adopted physical 
conditions specified in the 1" entry of Table 2 in WC98). 

To test our chemical kinetic code we have adopted rate 
constants only for the forward reaction^ (282 - 136/2+214- the 
network of WC98 requires the rate equations for 564 reactions) 
and we have computed the rate constants for all (282) backward 
reactions according to the method described above (hereafter 
case 1). As expected, in this case our code reproduces the WC98 
equilibrium concentrations. The results are marked as crosses 
in Fig.|6] The species (except atoms collected at the right side 
of the plot) are plotted in order of decreasing abundances. In 
addition, we have tested stability of our code by assuming that 
initial chemical composition of all species is equal to that at the 
thermodynamical equilibrium for the considered conditions. No 
changes in the species concentrations were seen after very long 
time (10** sec). Therefore, for further computations we used the 
equilibrium concentrations as an input to our code (for its faster 
convergence). Note that the results are the same when we start 
with a purely atomic gas. 

In the next test (hereafter case 2) we took all rate constants 
available in the WC98 network (i.e. 350=214+136 - see the 
previous subsection) and computed the rate constants only 
for the missing (214) backward reactions by using NIST 
thermochemical data (in fact, this approach is the same as 
described by WC98). The results are presented in Fig.|6] as 
circles. It is seen, that the present steady state differs from 
the thermodynamical equilibrium concentrations (case 1). More 
detailed comparison for all species, which belong to a given 
family (molecules containing element: H, C, N, O, Si, and S, 
respectively) are presented in six panels in Fig. as circles: 
going from the left to right and from the top to bottom for: 



* As a forward reaction, among reactions with rate constants given 
for both directions, we chose the reaction with more precise rate 
constant i.e. first when Arrhenius parameter /3 (see Eq. II il l is different 
from zero. If /? = for both reactions, we took reaction with smaller 
Arrhenius parameter In case of reactions without the backward rate 
constants we treat all of them as a forward reactions. 



H, C, N, O, Si, and S. For each panel we defined the separate 
ordinal number of species according to decreasing equilibrium 
concentrations with exception of atom(s) which are shown on 
the right side of each panel. 

The inconsistency in species abundances is rather high and 
only fifteen species differ less then 30 % (i.e. He, CO, H2, SiH, 
SiH3, S, HCO, Si, H, CS, SiHi, SiH4, HS, SiS, N2, and N). 
The highest differences are seen for oxygen-based molecules. 
Only formyl (HCO) and the most abundant among O-bearing 
molecules - carbon monoxide (CO) are in relatively good 
agreement with the corresponding equilibrium concentrations. 
Concentrations of other molecules (i.e. SiO, H2O, OH, H2CO, 
CO2, NO, SO, and O2) are far from equilibrium values. 

The WC98 network was created on the basis of the RATE95 
database (iMillar et alj|1997h . The discrepancies present in the 
above results encourage us to update the rate constants using 
the UDFA05 database (IWoodall et all I2OO6I) . We found that 
162 reactions among 564 have now the new values of the rate 
constants. We performed the same test as in case 2 (i.e. only 
missing backward rates were computed using Eq.[l4]l with the 
reaction network updated using the UDFA05 database. The 
results are presented as squares in Fig.|6] and the family splitting 
in this case is shown in Fig. |#lalso by squares. 

The agreement with the equilibrium state for oxygen-bearing 
molecules now is much better However, the discrepancies 
in abundances of carbon-bearing molecules are still present. 
Moreover, the differences with equilibrium results for nitrogen- 
bearing molecules are even higher 

The most likely source of inconsistency is the problem with 
the determination of reaction rate constants for the silicon- 
and sulfur-bearing species. WC98 estimated the rate constants 
for these reactions from similar reactions involving covalent 
species, i.e. oxygen for sulfur and carbon for silicon (see 
comments in Table 5 of WC98). When this group of reactions 
was excluded from the studied network, the abundances of 
carbon-bearing and oxygen-bearing molecules became more 
consistent with LTE results. However, the results for nitrogen- 
bearing species are still unsatisfactory and additional work on 
the chemical network extension is needed. 



4. Conclusions 

In this paper we have performed tests of our chemical 
(equilibrium and kinetic) codes with the aim to use them during 
self-consistent modelling of dynamics and chemistry in outflows 
from C-rich AGB stars. 

The LTE molecular distribution at the base of the wind 
is needed as the boundary condition to determine chemical 
composition in the circumstellar envelope of the AGB star. 
In such LTE calculations it is necessary to keep in mind 
the possible influence of the different sets of equilibrium 
dissociation constants on the species concentration. Analysis 
of the conducted test calculations shows that, fortunately, such 
influence for not abundant molecules is rather small. However, 
for more abundant species the influence is crucial and the used 
sets of equilibrium dissociation constants should be carefully 
selected and checked. In particular, our calculations have shown 
that the NIST database reproduce well the WC98 equilibrium 
concentrations, while agreement between WC98 and Tsuji is 
much worse. 

Steady state solution obtained with kinetic code for 
WC98 reaction network is different from the thermodynamical 
equilibrium solution. This would affect the full time-dependent 
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radiative-hydrodynamic computations of the inner wind. Note 
that CN and C2 molecules, important for the opacity 
computations in the atmosphere of carbon rich star, are 
underabundant by more than one order of magnitude in 
comparison to equilibrium concentration (see Fig|6]l. Our kinetic 
computations also show the strong overabundance of oxygen- 
bearing molecules (especially of SiO, H2O, and OH) in 
comparison to the LTE approach. It would be interesting to 
investigate how these overabu ndances influenced the final re sults 
of the in ner wind studies by IWillacv & Cherchnefll (Il998l) and 
ICherchne ff (2006). 

To make both, LTE and kinetic, models consistent we 
propose to replace all reaction rates in backward directioi£| by 
reaction rates computed from forward reactions with usage of 
thermochemical data. Then, in the limit of high density and 
temperature (i.e. when chemical timescales are much shorter 
than dynamical timescales) the kinetic steady state solution 
approach the local thermodynamic al e quilibrium. 

Observations dTsuii et al.l |1973|) and the astrochemical 



modefling (WC98. ICherchnefil l2006l) are showing that the 



base of inner wind in AGBs is a region of active chemistry. 
Therefore, chemical network, which may be used to investigate 
the composition of this region should include extended set of 
chemical reactions such as neutral-neutral involving metals and 
photochemistry induced by the stellar radiation field. This will 
be direction of our future studies. 
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